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Fast unsupervised Bayesian image 
segmentation with adaptive spatial 
regularisation 
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Abstract 

This paper presents a new Bayesian estimation technique for hidden Potts-Markov random 
fields with unknown regularisation parameters, with application to fast unsupervised A'-class image 
segmentation. The technique is derived by first removing the regularisation parameter from the 
Bayesian model by marginalisation, followed by a small-variance-asymptotic (SVA) analysis in 
which the spatial regularisation and the integer-constrained terms of the Potts model are decoupled. 

The evaluation of this SVA Bayesian estimator is then relaxed into a problem that can be computed 
efficiently by iteratively solving a convex total-variation denoising problem and a least-squares clus¬ 
tering (K-means) problem, both of which can be solved straightforwardly, even in high-dimensions, 
and with parallel computing techniques. This leads to a fast fully unsupervised Bayesian image 
segmentation methodology in which the strength of the spatial regularisation is adapted automatically 
to the observed image during the inference procedure, and that can be easily applied in large 2D 
and 3D scenarios or in applications requiring low computing times. Experimental results on real 
images, as well as extensive comparisons with state-of-the-art algorithms, confirm that the proposed 
methodology offer extremely fast convergence and produces accurate segmentation results, with the 
important additional advantage of self-adjusting regularisation parameters. 
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Image segmentation, Bayesian methods, spatial mixture models, Potts Markov random field, 
convex optimisation. 


I. Introduction 


Image segmentation is a eanonieal inverse problem whieh involves elassifying image pixels 
into elusters that are spatially eoherent and have well defined boundaries. It is widely aeeepted 
that this task ean be formulated as a statistieal inferenee problem and most state-of-the- 
art image segmentation methods eompute solutions by performing statistieal inferenee (e.g., 
eomputing penalized maximum likelihood or maximum-a-posteriori estimates). In this paper 
we foeus on new Bayesian eomputation methodology for hidden Potts-Markov random fields 
(MRFs) 0. a powerful elass of statistieal models that is widely used in Bayesian image 
segmentation methods (see 0-0 for reeent examples in hyperspeetral, non destruetive 
testing, ultrasound, and fMRI imaging). 

Despite the wide range of applieations, performing inferenee on hidden Potts MRFs re¬ 
mains a eomputationally ehallenging problem. In partieular, eomputing the maximum-a- 
posteriori (MAP) estimator for these models is generally NP-hard, and thus most image 
proeessing methods eompute approximate estimators. This has driven the development of 
effleient approximate inferenee algorithms, partieularly over the last deeade. The eurrent 
predominant approaehes for approximate inferenee on MRFs are based on eonvex models 
and eonvex approximations that ean be solved effieiently by eonvex optimisation @-[[8|, and 
on approximate estimators eomputed with graph-eut Q, and message passing algorithms 
l[TT|-[[Tg. In a similar fashion, modern algorithms to solve aetive eontour models, the other 
main elass of models for image segmentation, are also prineipally based on eonvex relaxations 


and eonvex optimisation [ 14|, [ 151 and on Riemannian steepest deseent optimisation sehemes 


An important limitation of these eomputationally effleient approaehes is that they are 
supervised, in the sense that require praetitioners to speeify the value of the regularisation 
parameter of the Potts MRF. However, it is well known that appropriate values for regularisa¬ 
tion parameters ean be highly image dependent and sometimes diffleult to seleet a priori, thus 
requiring praetitioners to set parameter values heuristieally or by visual eross-validation. The 
Bayesian framework offers a range of strategies to eireumvent this problem and to design 
unsupervised image segmentation inferenee proeedures that self-adjust their regularisation 
parameters. Unfortunately, the eomputations involved in these inferenees are beyond the 
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scope of existing fast approximate inferenee algorithms. As a eonsequenee, unsupervised 
image segmentation methods have to use more eomputationally intensive strategies sueh 


as Monte Carlo approximations [|^|, [211, variational Bayes approximations [|S|, and EM 
algorithms based on mean-field like approximations @1. 

In this paper we propose a highly effieient Bayesian eomputation approaeh speeifieally 
designed for performing approximate inferenee on hidden Potts-Markov random fields with 
unknown regularisation parameters, with applieation to fast unsupervised A"-elass image 
segmentation. A main originality of our development is to use a small-varianee-asymptotie 
(SVA) analysis to design an approximate MAP estimator in whieh the spatial regularisation 
and the integer-eonstrained terms of the Potts model are deeoupled. The evaluation of this 
SVA Bayesian estimator ean then be relaxed into a problem that ean be eomputed effieiently 
by iteratively solving a eonvex total-variation denoising problem and a least-squares elustering 
(K-means) problem, both of whieh ean be solved straightforwardly, even in high-dimensions, 
and with parallel eomputing teehniques. 


Small-varianee asymptoties estimators were introdueed in [251 as a eomputationally effi¬ 
eient framework for performing inferenee in Diriehlet proeess mixture models and have been 
reeently applied to other important maehine learning elassifieation models sueh as the Beta 
proeess and sequential hidden Markov models [|26|, as well as to the problem of configuration 


alignment and matehing [27|. Here we exploit these same teehniques for the hidden Potts 
MRF to develop an aeeurate and eomputationally effieient image segmentation methodology 
for the fully unsupervised ease of unknown elass statistieal parameters (e.g., elass means) 
and unknown Potts regularisation parameter. 

The paper is organised as follows: in Seetion II we present a brief baekground to Bayesian 
image segmentation using the Potts MRF. This then followed by a detailed development 
of our proposed methodology. In Seetion IV the methodology is applied to some standard 
example images and eompared to other image segmentation approaehes from the state of the 
art. Finally some brief eonelusions are drawn in Seetion V. 


II. Background 

We begin by reealling the standard Bayesian model used in image segmentation problems, 
whieh is based on a finite mixture model and a hidden Potts-Markov random field with known 
regularisation parameter ( 3 . For simplieity we foeus on univariate Gaussian mixture models. 
However, the results presented hereafter ean be generalised to all exponential-family mixture 
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models (e.g., mixtures of multivariate Gaussian, Rayleigh, Poisson, Gamma, Binomial, ete.) 
by following the approaeh deseribed in 

Let r/n G M denote the nth observation (i.e. pixel or voxel) in a lexieographieal veetorized 
image y = {yi, ^ ■ We assume that y is made up by K regions {Ci, ..., Ck} 

sueh that the observations in the fcth elass are distributed aeeording to the following eondi- 
tional marginal observation model 

r/n|n G Cfc ~ A/'(/ifc,cr^), (1) 

where G M represents the mean intensity of elass Ck- For identifiability we assume that 
^^k 7 ^ l^j for all k ^ j. 

To perform segmentation, a label veetor 2 = {zi, ..., zn) is introdueed to map or elassify 
observations y to elasses Ci, ... ,Ck Zn = k if and only if n G Ck). Assuming that 
observations are eonditionally independent given 2 and given the parameter veetor y, = 
(/xi, ..., yx), the likelihood of y ean be expressed as follows 

K 

fiy\z,y) = Yl Ylp^ivnlyk^cr"^), (2) 

k=l n£Sk 

with Sk = {n : Zn = k}. A Bayesian model for image segmentation is then defined by 
speeifying the prior distribution of the unknown parameter veetor {z, y). The prior for 2 is 
the homogenous iL-state Potts MRF p9| 

fY\l3) = ^ 7 ^ exp [l3H{z)l (3) 

with regularisation hyper-parameter /? G M"*", Hamiltonian 

N 

H{z) = '^ 

n=l n'eV(n) 

where 5(-) is the Kroneeker funetion and V(?x) is the index set of the neighbors of the rxth 
voxel (most methods use the 1st order neighbourhoods depieted in Fig. [^, and normalising 
eonstant (or partition funetion) 

= X^exp[/iiT(2)]. (5) 

Z 

Notiee that the Potts prior Q is defined eonditionally to a given value of jS. Most image 
segmentation methods based on this prior are supervised; i.e., assume that the value of f3 is 
known and speeified a priori by the praetitioner. Alternatively, unsupervised methods eonsider 
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that /3 is unknown and seek to adjust its value automatieally during the image segmentation 
proeedure (this point is explained in detail in Seetion ??). 

In a similar fashion, the elass means are eonsidered prior independent and assigned Gaus¬ 
sian priors /i^ ~ A^(0,p^) with fixed varianee p^, 

K 

= YlPM{f^k\0,p‘^)- ( 6 ) 

k=l 

(to simplify notation the dependenee of distributions on the fixed quantity is omitted). 

Then, using Bayes theorem and taking into aeeount the eonditional independenee strueture 
of the model (see Fig. [^, the joint posterior distribution of ( 2 ,/i) given y and /3 ean be 
expressed as follows 

cx f{y\z,y)f{z\l3)f{y), (7) 

where oc denotes proportionality up to a normalising eonstant that ean be retrieved by setting 
J / (z, /3) dzd^i = 1. The graphieal strueture of this Bayesian model is summarised in 

Fig. below. Notiee the Markovian strueture of 2 and that observations are eonditionally 
independent given the model parameters 2 , pi and cr^. 



Fig. 1. [Left:] Directed acyclic graph of the standard Bayesian model for image segmentation (parameters with fixed values 
are represented using black boxes). [Right] Local hierarchical representation of the hidden Potts MRF and the observed 
image for 4 neighbouring pixels. 


Finally, given the Bayesian model 0, a segmentation of y is typieally obtained by 
eomputing the MAP estimator 

zi,/ii = argmax / (2,^|f/,/3), (8) 

Z,fJL 
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Fig. 2. 4-pixel (left) and 6-voxel (right) neighborhood structures. The pixel/voxels considered appears as a void red circle 
whereas its neighbors are depicted in full black and blue. 


which can also be obtained by solving the equivalent optimisation problem 

= argmin-log/. (9) 

Unfortunately these optimisation problems are known to be NP-hard due to the combinatorial 
nature of the Potts Hamiltonian H{z) defined in (Q. As mentioned previously, modern 
image segmentation methods based on (|^ typically address this issue by using approximate 
(local) integer optimisation algorithms (e.g., graph-cut, message passing) (ig-iiTg, and more 
recently with convex relaxations of the Potts model (see for instance [[^, [(7|). 


III. Proposed method 

This section presents a highly computationally efficient approach for performing approx¬ 
imate inference on z when the value of the regularisation parameter /) is unknown. The 
approach is based on a small-variance asymptotics (SVA) analysis combined with a convex 
relaxation and a pseudo-likelihood approximation of the Potts MRF. Our development has 
three main steps. In the first step we adopt a hierarchical Bayesian approach to remove /) 
from the model by marginalisation; because marginalising w.r.t. fi requires knowledge of the 
intractable Potts partition function (|^ we use a pseudo-likelihood approximation. However, 
performing inference with the resulting marginalised model is still NP-hard. In the second part 
of our development we address this difficulty by using auxiliary variables and an SVA analysis 
to decouple the spatial regularisation and the integer-constrained terms of the Potts model. 
The evaluation of the resulting SVA Bayesian estimator is then relaxed into a problem that 
can be computed efficiently by iteratively solving a convex total-variation denoising problem 
and a least-squares clustering problem, both of which can be solved straightforwardly, even in 


high-dimensions, with parallel implementations of Chambolle’s optimisation algorithm [30| 
and of K-means gg. 
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A. Marginalisation of the regularisation parameter {3 

Following a hierarchical Bayesian approach, we address the fact that the value of (3 is 
unknown by modelling it as an additional random variable of the Bayesian model. Precisely, 
we assign [3 a prior distribution /(/3) and define an augmented model that includes f3 within its 
unknown parameter vector. By using Bayes’ theorem we obtain the joint posterior distribution 

f {x,z,p,l3\y) oc f{y\x)f{x\z,p)f{p)f{z\l3)f{l3) (10) 

which includes [3 as an unknown variable. The rationale for replacing the fixed regularisation 
parameter /? of 0 by a random variable with prior f{(3) is that it is often possible to specify 
this prior distribution such that the amount of regularisation enforced by the Potts MRF 
is driven by data and the impact of f{(3) on the inferences is minimal. At the same time, 
experienced practitioners with knowledge of good values of /3 can specify f{(3) to exploit 
their prior beliefs. In this paper we use a gamma (hyper-)prior distribution 

/(/3) = 7 “/^“-'exp (-7/?)lM+(/3)/r(a) 


because it has favourable analytical tractability properties that will be useful for our devel¬ 
opment (appropriate values for the fixed parameters a and 7 will be derived later through a 
small-variance asymptotics analysis). 

Moreover, in order to marginalise (3 from the model we notice that (3 is conditionally 
independent of y given z; to be precise, that / {x, z, p,, I3\y) = f {f3\z)f{x, z, p\y). There¬ 
fore, integrating f {x, z, p, (3\y) with respect to /3 is equivalent to redefining the posterior 
distribution ( [T^ with the marginal prior /(z) = f{z, /3)df3. Evaluating this marginal 
prior exactly is not possible because it requires computing the normalising constant of the 


Potts model C{(3) defined in 0, which is a reputedly intractable problem [20|. To obtain 
an analytically tractable approximation for this marginal prior we adopt a pseudo-likelihood 


approach |32| and use the approximation C{(3) oc (3 leading to 

f{z)= [ /(2,/3)d/3 
Jr+ 

oc [ /3^exp(/3i7(z))/3“-iexp(-7/?)d/? 

Jr+ 

and to the following (marginal) posterior distribution 

K 

f{x,z,p\y) oc /(^)( 7 -Ff(z))“^"+^^ JJ JJ {yn\pk,(r^) , 

k=l n&Sk 

that does not depend on the regularisation parameter f3. 


( 11 ) 


( 12 ) 
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B. Small-variance approximation 


The next step of our development is to eonduet a small-varianee asymptoties analysis on 
( fT^ and derive the asymptotic MAP estimator of a;, 2 , p. We begin by introducing a carefully 
selected auxiliary vector x such that y and {z, p) are conditionally independent given x, 
and that the posterior / {x, z, p\y) has the same maximisers as (|7]) (after projection on the 
space of (z, p)). More precisely, we define a random vector x G with degenerate prior 

K 

f{x\z,p) = W_ Yl6{xn-Pk), (13) 

k=l n£Sk 

and express the likelihood of y given x, z and p as 

N 

f{y\x,z,p) = f{y\x) = Y[pM(,yn\xn,(J^)- 

n=l 

The prior distributions for 2 and p remain as defined above. The joint posterior distribution 
of X, 2 , p is given by 


f {x,z,p,l3\y) oc f{y\x)f{x\z,p)f{z\/3)f{p) 


oc 


K 


n n PjviVnlXn, 0-'^)S{Xn - Pk) 

k=l n&Sk 


[7 - . 


(14) 


Notice that from an inferential viewpoint ( pA] ) is equivalent to ( [T^ , in the sense that marginal¬ 
ising X in ( [T4l ) results in ( [T^ . 

Moreover, we define H*{z) as the “complement” of the Hamiltonian H{z) in the sense 
that for any 26 [ 1 , • • •, K]^ 

H{z) + H*{z) = N\V\, 


where |V| denotes the cardinality of the neighbourhood structure V. For the Potts MRF this 
complement is given by 

N 

(15) 

n=l n'eV(n) 

Replacing H{z) = N\V\—H*{z) in ( [T4| ) we obtain 

f {x,z,p,P\y) Yl PJ^iyn\xn,(^^)S{xn - PkU f{p)[H*{z) + {-f - (16) 

Furthermore, noting that H*{z) only measures if neighbour labels are identical or not, 
regardless of their values, it is easy to check that the posterior ( [T^ remains unchanged if we 
substitute H*{z) with H*{x) 


K 

f {x,Z,p,l3\y) CX f{p) [H*{x) + (7 - JJ JJ PN{yn\Xn,(T‘^)5{Xn - Pk)- (17) 

fc=l n&Sic 
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Finally, we make the observation that for 1st order neighbourhoods (see Fig. we have 
H*{x) = 2||Va;||o, where ||Va;||o= ||V/ia;||o+||V^a;||o denotes the Iq norm of the horizontal 
and vertieal eomponents of the 1 st order diserete gradient of x, and therefore 

K 

f {x,z,ii,P\y) oc f{y) [||Va;||o+(7 - iV|V|)/2]“^“+^^ JJ JJ PM{yn\xn,(y‘^)5{xn - yu)- (18) 

k=l n&Sk 

The graphieal strueture of this equivalent hierarehieal Bayesian model is summarised in 
Fig. I^below. Notiee that in this model x separates y and cx^ from the other model parameters, 
that the regularisation parameter (5 has been marginalised, that the MRF is now enforeing 
spatial smoothness on x not z, and that the elements of z are prior independent. 



Fig. 3. [Left:] Directed acyclic graph of the proposed Bayesian model, augmented by the auxiliary variable x decoupling 
/I and 2 from y, and with marginalisation of the regularisation parameter j3 (parameters with fixed values are represented 
using solid black boxes, marginalised variables appear in dashed boxes). [Right] Local representation of three layers of the 
model for 4 neighbouring pixels. 


We are now ready to eonduet a small-varianee asymptotics analysis on ([181) ^nd derive the 


asymptotic MAP estimator of x,z,iJ,, which is defined for our model as [25| 


argmin lim —a'^ log f {x, z, y,\y) . 

x,z,y a-'^^o 
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First, we use the faet that 5(s) = lim^ 2 _^oPAr(s| 0 ,r^) to express ( [T8] ) as follows 
f {x,z,^i\y,P) 


K 


oc 


, n n PAfiVnlXn, (j‘^)PAr{Xn\Pk, ' 


t2-).0 


^^=1 n£Sk 


x/(/>) |l|Va:||„+(7-lV|V|)/2] 


a^N 


(19) 


K 


OC 


.‘i“„ n n 


T-^^O 


exp 


(Xn - Vnf (Xn “ Pkf 


2a2 


\k=l riGSk 

x/(^)[||Vir||o+(7-iV|V|)/2]-(“+^) 


2r2 


Then, in a manner akin to Broderiek et al. [251, we allow the model’s hyper parameters to 
seale with cr^ in order to preserve the balanee between the prior and the likelihood and avoid 
a trivial limit. More precisely, we set a = N/a'^ and assume that vanishes at the same 
speed as Then, the limit of —log / {x, z, n\y) as —)■ 0 is given by 


K 


lim -CT^ log/(ai,Z,^|t/) - Vnf + ^{Xn - Pkf 


O-2_s.0 


k=l n£Sk 

+ iVlog(||Vir||o+(7-iV|V|)/2), 

and the MAP asymptotic estimators of x,z, fj, by 


( 20 ) 


K 


argmin ^ ^ hxn - Vnf + ^(ain - Pkf + log(||Vai||o+l), 

k=i nes, ^ 


( 21 ) 


where we have set 7 = 2 + A^|V| such that the penalty log [||Vai||o+(7 — A^|V|)/2] > 0. 


C. Convex relaxation and optimisation 


Computing the estimator pTj ) is still NP-hard due to log(| |Vai| |o+l). To address this 
difficulty we use a convex relaxation of ||Vai||o and exploit the concavity of the logarithmic 
function. Precisely, we replace 11 Vai| |o by the convex approximation TV(ai) = 11 Vai| |i_ 2 , (i.e.. 


the isotropic total-variation pseudo-norm of x [331), and obtain the following optimisation 
problem 


K 


argmin^ Vnf+ \{.Xn-Pkf + N\og{TV{x)+ 1), (22) 




k=l 


which can be very efficiently computed by iterative minimisation w.r.t. x, z and p. The 
minimisation of ( [22] ) w.r.t. 2 (with x and p fixed) is a trivial separable integer problem that 
can be formulated as N independent (pixel-wise) minimisation problems over 1,..., iC (these 
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unidimensional integer problems ean be solved by simply ehecking the value = 1,..., 
that minimises ( |22l ) for eaeh pixel n = 1,..., iV). Similarly, the minimisation with respeet 
to ^ is a trivial quadratie least squares fitting problem with analytic solution (i.e., by setting 
^J^k = \^\ ^ ^ - ■ ■ iK, where \Sk\ denotes the cardinality of Sk)- Also 

note that iteratively minimising ( |22l ) with respect to z and with fixed x, is equivalent 
to solving a least squares clustering problem with the popular K-means algorithm [[3T|. 


Moreover, the minimisation of ( [22| ) w.r.t. x (with z and fixed) is achieved by solving 
the non-convex optimisation problem 


K 


argmin ^ - Vnf+ ^{xn - + Nlog[TV{x) + 1], 


(23) 


k=l n^Sk 


which was studied in detail in [34|. Essentially, given some initial condition E ( |2^ 
can be efficiently minimised by majorisation-minimisation (MM) by iteratively solving the 
following sequence of trivial convex problems. 


V 


^ 1 1 

(^+1) = argmin ^ ^ -{xn - |/n)^ + -{xn - /J-kf + XeTV{x), 

k=l n£Sk 


(24) 


with Xf = 


N 


ri/[uW] + 1’ 

in which plays the role of a regularisation parameter, and where we have used the majorant 


[341 




(25) 


> log {TV{v^^^) + 1) . 

Notice that each step of ( |24| ) is equivalent to a trivial convex total-variation denoising 
problem that can be very efficiently solved, even in high-dimensional scenarios, by using 
modern convex optimisation techniques (in this paper we used a parallel implementation of 


Chambolle’s algorithm [30|). 


The proposed unsupervised segmentation algorithm based on ( |22| ) is summarised in Algo, 
below. We note at this point that because the overall minimisation problem is not convex 


the solution obtained by iterative minimisation of ( |22| ) might depend on the initial values of 
X, z, fx. In all our experiments we have used the initialisation = 2y, z = [1,..., 1]^, 
= [0,..., 0]^ that produced good estimation results. 
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Algorithm 1 Unsupervised Bayesian segmentation algorithm 
1: Input: Image y, number of maximum outer iterations T and inner iterations L, toleranee 

level e. 

2 : Initialise = 2y, z = [1,..., 1]^, y, = [0,..., 0]^. 

3: for t = 1 : T do 
4: Set 

5: for f = 0 : L do 

6: Set = N/{TV[v^^^] + 1}. 

7: Compute using ( |24| ), with fixed z = and y = using Chambolle’s 

algorithm [ [30| . 

8: if {N/{TV[v^^+^^] + 1} - A) > eA then 

9: Set £ = £+!. 

10 : else 

11 : Exit to line 14. 

12 : end if 

13: end for 

14: Set 

15: Compute and y^*'^ by least-squares elustering of x^*'> using the K-means algorithm 

m- 

16 : if z^*) 7 ^ then 

17: Set t = t + 1. 

18 : else 

19: Exit to line 22. 

20 : end if 

21 : end for 

22 : Output: Segmentation y^^\ A = N/{TV[x^^'>] + 1). 


IV. Experimental Results and Observations 


In this seetion we demonstrate empirieally the proposed Bayesian image segmentation 
methodology with a series of experiments and eomparisons with state-of-the-art algorithms. 
To asses the aeeuraey of our method we eompare the results with the estimations produeed by 


the Markov ehain Monte Carlo algorithm [20|, whieh estimates the marginal posterior of the 
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segmentation labels f {z\y) with very high aeeuraey. We also report eomparisons with four 
supervised fast image segmentation teehniques that we haven ehosen to represent different 
effieient algorithmie approaehes to image segmentation (e.g. MRF energy minimisation solved 
by graph-eut, aetive eontour solved by Riemannian gradient deseent, and two convex models 
solved by convex optimisation). The specific methods used in the comparison are as follows: 

• The two-stage smoothing-followed-by-thresholding algorithm (TSA) p3| , which is closely 
related to a semi-supervised instance of Algo. with a single iteration (TV-denoising 
followed by K-means), and with a fixed regularisation parameter A specified by the 
practitioner. 

• Hidden Potts MRF segmentation Q with fixed (3, solved by graph max-flow/min-cut 


approximation [351. 


Chan-Vese active contour by natural gradient descent [ 16| (to our knowledge this method 
is currently the fastest approach for solving active contour models). 


• The fast global minimisation algorithm (FGMA) [14| for active contour models. In a 
similar fashion to our method, this algorithm also involves a model with a TV convex 
relaxation that is solved by convex optimisation. 

We emphasise that, unlike the proposed method, all these efficient approaches are supervised, 
i.e., they require the specification of a regularisation parameters. In the experiments reported 
hereafter we have tuned and adjusted the parameters of each algorithm to each image by 
use of visual cross-validation to ensure we produce the best results for each method on each 
image. 

To guarantee that the comparisons are fair we have applied the six algorithms considered 
in this paper to three images with very different compositions: the Lungs and Bacteria 
images from the supplementary material of [ |T4| , and one slice of a 3D in-vivo MRI image of 
a human brain composed by biological tissues (white matter and grey matter) with complex 
shapes and textures, making the segmentation problem challenging. The three test images 
are depicted in Figure These images have been selected as they are composed of different 
types and numbers of objects; objects which have different shapes, (regular and irregular); 
and a range of potential segmentation solutions. All experiments have been conducted using 
a MATLAB implementation of Algo. [I] with parameters T = 50, L = 25, e = 10“^, and 
computed on an Intel i7 quad-core workstation running MATLAB 2014a. With regards to the 
algorithms used for comparison, when possible we have used MATLAB codes made available 
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by the respective authors. It should be noted that these are mainly MATLAB scripts, however 


the graph-cut method is written in C-i-i-, ( the [361 implementation was used here), so it has 
a slight advantage in terms of computational performance. 

We emphasise at this point that we do not seek to explicitly compare the accuracy of 
the methods because: 1) there is no objective ground truth; 2) the ’’correct” segmentation 
is often both subjective and application-specific; and 3) the segmentations can often be 
marginally improved by fine tuning the regularisation parameters. What our experiments 
seek to demonstrate is that our method performs similarly to the most efficient deterministic 
approaches of the state-of-the-art, both in terms of segmentation results and computing 
speed, with the fundamental advantage that it does not require specification of the value 
of regularisation parameters (i.e., it is fully unsupervised). 

Figures and [^respectively show the segmentation results obtained for the Lungs, 
Bacteria and Brain test images with each method. The segmentations of the Lungs 
and Bacteria images have been computed using K = 2 classes to enable comparison 


with the natural gradient method [ 16| and FGMA [ 14 j (these methods are based on an active 
contour model that only supports binary segmentations), whereas the Brain image has been 
computed using K = 3 classes to produce a clear segmentation of the grey matter and the 
white matter. The computing times associated with these experiments are reported in Table 
|I| Observe that all six methods produced similar segmentation results that are in good visual 
agreement with each other. In particular, we observe that the proposed method successfully 
determined the appropriate level of regularisation for each image and produced segmentations 


that are very similar to the results obtained with the supervised methods graph-cut [351 
and TSA [15[, and with the unsupervised MCMC algorithm [ [T^ that in a sense represents 
a benchmark for these approximate inference methods. Moreover, Table |I| shows that the 
proposed method was only 2 or 3 times slower than state-of-the-art supervised approaches, 
which is an excellent performance for a fully unsupervised method. This additional computing 


time is mainly due to the additional computations related to the non-convex program ( [23] ); 
however, we emphasise that this algorithm has the property of adapting automatically the 
level of regularisation to the image, and that the computing times reported in Table |I| do not 
take into account the time involved in running the supervised algorithms repeatedly to adjust 
their regularisation parameters. 
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TABLE I 

Computing times (seconds) eor the lungs, bacteria and brain images displayed in Figs.[^ Figs.[^and 

Figs.[71 



Bacteria 

Bacteria 

Brain 

Proposed 

0.65 

0.80 

0.23 

TSA 


0.20 

0.21 

0.17 

Graph-Cut 

351 

0.30 

0.30 

0.21 

Natural gradient 

16 

0.20 

0.18 

n/a 

FGMA 1 

14| 

0.32 

0.47 

n/a 

MCMC 

16 


900 

1150 

533 


V. Conclusions 

We have presented a new fully unsupervised approaeh for eomputationally effieient image 
segmentation. The approaeh is based on a new approximate Bayesian estimator for hidden 
Potts-Markov random fields with unknown regularisation parameter (3. The estimator is based 
on a small-varianee-asymptotie analysis of an augmented Bayesian model and a eonvex 
relaxation eombined with majorisation-minimisation teehnique. This estimator ean be very 
effieiently eomputed by using an alternating direetion seheme based on a eonvex total- 
variation denoising step and a least-squares (K-means) elustering step, both of whieh ean be 
eomputed straightforwardly, even in large 2D and 3D seenarios, and with parallel eomputing 
teehniques. Experimental results on real images, as well as extensive eomparisons with 
state-of-the-art algorithms showed that the resulting new image segmentation methodology 
performs similarly in terms of segmentation results and of eomputing times as the most 
effieient supervised image segmentation methods, with the important additional advantage of 
self-adjusting regularisation parameters. A detailed analysis of the theoretieal properties of 
small-varianee-asymptotie estimators in general, and in partieular of the methods deseribed 
in this paper, is eurrently under investigation. Potential future researeh topies inelude the 
extension of these methods to non-Gaussian statistieal models from the exponential family 
and their applieation to ultrasound and PET image segmentation, extensions to models with 
unknown number of elasses K, and eomparisons with other Bayesian segmentation methods 
based on alternative hidden MRE models that ean also be solved by eonvex optimisation, 
sueh as d- 
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(b) Bacteria 
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(c) Brain 

Fig. 4. The Lungs (336 x 336 pixels), Bacteria (380 x 380 pixels), and Brain (256 x 256 pixels) images 
used in the experiments. 
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(a) Proposed 


(b) MCMC 1201 




(C) TSA 


(f) Graph-Cut |35| 



(e) Natural grad, 


(f) FGMA 1141 


Fig. 5. Comparison with the state-of-the-art methods 1161, p5) , p4) , and GD using the lung image (336 x 336 


pixels) from the supplementary material of 
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Fig. 6. Comparison of the supervised and unsupervised methods with the state of the algorithm | [T6) , p5) , | [T4) 

and p3) using the bacteria image (380 x 380 pixels) from the supplementary material of | [T4) . 
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(C) TSA ^ 

Fig. 7. Segmentation of a brain MRI image (256 x 256 pixels). 


(d) Graph-Cut |35| 
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